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Abstract The growth of evolutionary computing (EC) methods in the exploration 
of complex potential energy landscapes of atomic and molecular clusters, as well 
as crystals over the last decade or so is reviewed. The trend of growth indicates 
that pure as well as hybrid evolutionary computing techniques in conjunction of 
DFT has been emerging as a powerful tool, although work on molecular clusters has 
been rather limited so far. Some attempts to solve the atomic/molecular Schrodinger 
Equation (SE) directly by genetic algorithms (GA) are available in literature. At 
the Born-Oppenheimer level of approximation GA-density methods appear to be a 
viable tool which could be more extensively explored in the coming years, specially 
in the context of designing molecules and materials with targeted properties. 
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1 Introduction 

An N-atom system can be arranged in many different ways to form chemically 
meaningful structures. These structures represent minima on the 3N dimensional 
potential energy surface (PES) on which the nucleii move. One of these minima 
is the global one and search for the global minimum is a computationally chal¬ 
lenging job. Optimization strategies play a very important role in theoretically elu- 
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cidating the minimum energy structures supported on a potential energy surface. 
The ideas of PES and minimum energy structures are rooted in the Born Oppen- 
heimer (BO) theory of molecules. BO theory starts by separating the slow nuclear 
motion from the much faster electronic motion which allows the electronic charge 
distribution in a molecule to adjust itself practically instantaneously as the nuclei! 
move. The separation results in the creation of the so called adiabatic potential en¬ 
ergy surface En{xi,X 2 ,‘ ■ ■ ,X 2 n) where £„s are the eigenvalues of a molecular elec¬ 
tronic Hamiltonian H{r,xi,X 2 , - ■ ■ ,X 2 n) in which the nuclear displacement coordi¬ 
nates (xi ,X 2 , • • • ,X 2 ,n) appear as parameters. The PES can have a number of station¬ 
ary points. At all such points 

dE„ 

^ =0, for/= 1,2,---,3A (1) 

OXi 

If in addition, the eigenvalues of the matrix of second derivatives of £„ with re¬ 
spect to 3A nuclear displacement coordinates (x,) are > 0, the stationary point is 
identified with a minimum energy structure. Out of 3A eigenvalues 3 represent¬ 
ing transnational and 3 representing rotational motion of the center of mass are 
discarded as motion along these coordinates do not define any new chemical struc¬ 
tures. It is immediately seen that the problem of locating minimum energy struc¬ 
tures on an N-atom PES is a problem of minimizing a 3N dimensional function 
En{Xi,X 2 , ■ ■ ■ ,X 3 n)- The function itself may be provided on the fly by solutions of 
the appropriate molecular Schrodinger equation or may be generated by superposing 
calibrated pair or higher body potentials. 

In addition to the minima on the PES, chemists are also interested in locating first 
order saddle points and constructing minimum energy path (the reaction path) that 
connects the saddle point to two neighboring minima. At the first order saddle 

dEn 

^ =0, for/= 1,2,---,3A (2) 

OXi 

and all but one of the eigenvalues (excluding the six) of the 2nd derivative matrix 
are positive, the special one being negative. 

The present review focuses attention on the ways of exploring PES of molecules, 
clusters, crystals, etc. by invoking techniques of evolutionary computing for locating 
minima. In addition, we also review attempts to solve the atomic or molecular elec¬ 
tronic Schrodinger equation (MESE) by similar techniques. In the case of molecular 
Schrodinger equation, the objective may be either to construct the PES by solving 
the MESE at a series of configuration and locating minima on the surface or to look 
for the deepest minimum in one go. 
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2 The Optimization Problem 


Any problem can in principle, be reduced to an equivalent optimization problem. 
The domain of the function to be optimized is called the search space (S). Goodness 
of any solution in S is measured with a mathematical function, called the objec¬ 
tive function or fitness function (O : S —> R). A fitness landscape (mapping from a 
configuration space into the real numbers) may be considered as a triple (S,0,D), 
where D is a metric defined on S. 

In general, optimization problems involve setting a vector X of free parameters of 
a system in order to optimize (maximize or minimize) some objective function 0(X) 
subject to the satisfaction of inequality constraints gi:(X), equality constraints hkO^), 
as well as upper and lower bounds on the range of allowable parameter values. 
Any constrained nonlinear optimization problem that deals with the search for a 
minimum of a nonlinear function 0(X) of m variables can be formulated as follows: 

minO(X), X = {x\, X 2 , x-i, ••• , Xmf € S (3) 

rgi(X)=0,k=l,2,...,n 

subject lo < hk{X) <0, k = (4) 

[ (x;)^ < Xi < {Xif 

Though equation Q is defined as minimization problem, it can be equivalently 
posed as a maximization problem by a simple modification of the objective func¬ 
tion, for example, by — 0(X), constant+o(x) several other means. There are 

no specific conditions attached here to the variable type and the function character¬ 
istics may have multimodal fitness landscapes and significant levels of parameter 
interactions. Other conditions which are commonly addressed in optimization re¬ 
search but not specifically addressed in the present review include dynamic objec¬ 
tive functions and multiple conflicting objectives. Let X S S be a feasible solution to 
the m—dimensional constrained nonlinear problem. X* be the optimum vector that 
solves equation The point X* S S is a local optimal solution to this problem, if 
there exists an e > 0 such that 


I 0(X) > 0(X*), VX G S : II X — X* II < e for local minimum 
|0(X) < 0(X*), VX G S : II X — X* II < e for local maximum 

X* G S is a global optimal solution if 

f 0(X) > 0(X*), VX G S for global minimum 
|0(X) < 0(X*), VX G S for global maximum 

The uniqueness of the global optimality is defined by 

f 0(X) > 0(X*), VX G S for global minimum 
[ 0(X) < 0(X*), VX G S for global maximum 


(5) 


(6) 


(7) 
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Global optimization methods differ from local optimization methods in that they 
attempt to find not just any local optimum, but the smallest (largest) local optimum 
in the search space S. Global optimization is a difficult problem since no general 
criterion exists for determining whether the global optimum has been reached. 


3 Types of Optimization Algorithms 

The goal of any optimization method is to design mathematical and computational 
infrastructure in order to locate the extremum (minimum or maximum) of a func¬ 
tion (or functional). Automated optimization algorithm designing (i.e. the algorithm 
can learn adaptively to tune all its parameters by its own) can be very costly and 
not always straightforward. Expertise is needed to tackle with the parameter sensi¬ 
tivity, efficiency, robustness and solvability of the algorithms. Even mathematically 
equivalent formulations often differ substantially in efficiency and solvability. So it 
requires careful thoughts along with mathematical details while designing these op¬ 
timization methods. Based on the use of random numbers, optimizers can be divided 
into two major categories, deterministic and stochastic. Deterministic algorithms are 
in general unidirectional (there exists at most one way to proceed, otherwise, the al¬ 
gorithm gets terminated) and do not use random numbers in any step of execution. 
State Space Search, Branch and Bound, Gradient based methods (such as Newton’s 
method, Gauss-Newton method, Steepest-Descent method, Levenberg-Marquardt 
algorithms). Conjugate - Direction methods (such as Conjugate-Gradient method, 
Eletcher-Reeves method, Powell’s method, Partan method), Quasi-Newton meth¬ 
ods (such as Davidon-Eletcher-Powell method, Broyden-Eletcher-Goldfarb-Shanno 
method, Hoshino method) are some examples of deterministic optimization methods 
m. The gradient-based methods are the most used classical techniques for solving 
optimization problems. However, these methods can only be applied to the objective 
functions that are continuous and differentiable. Some of them, such as the Newton 
method, even require the knowledge of the second derivatives of objective functions 
as well. These techniques are also limited to the optimization problems having a 
sole optimum, as they are gradient-based. However, most of the objective functions 
in real-world problems are not differentiable and they have a large number of local 
optima, which makes these methods inapplicable. The unavailability of analytical 
gradients can be tackled by using numerical gradients. But this is only feasible for 
relatively low dimensional problems. So gradient-free optimization procedures are 
often much sought after techniques in fairly large scale optimization. 

Stochastic algorithms, on the other hand, incorporate the concept of probability, 
employs at least one instruction or at least one operation that makes use of ran¬ 
dom numbers and do not use the gradient or Hessian matrix. The function to be 
optimized need not even be continuous or differentiable. If the relation between a 
candidate solution and its “fitness” are not so obvious or too complicated, or the 
dimensionality of the search space is very high, it becomes harder to solve a prob¬ 
lem deterministically. Apart from these reasons there can be many other pitfalls and 
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booby traps EEllllEllTlEllSl (such as conflicting objectives, heavily constrained 
fitness function, oversimplification, overfitting, non-separability, scalability, rugged 
and deceptive fitness landscape, neutrality, epistasis etc.) that make the optimiza¬ 
tion problems difficult to handle and can lead an optimizer to a suboptimal region 
in the search space. Simulated annealing, Monte-Carlo sampling. Stochastic tunnel¬ 
ing, Parallel tempering. Stochastic Hill Climbing, PSO, GA, Evolution Strategies, 
Memetic Algorithms, Differential Evolution are some examples of stochastic algo¬ 
rithms iiii [ini [ni [II [i3] El |3] a |5] [Ml i6i H1 [la . Our focus in this review has been 
mainly on pure evolutionary computing techniques and their hybrids in which power 
of several methods are combined to produce a superior algorithm for global search. 
Averaged over the set of all problems, any pair of optimizers perform equivalently, 
which essentially implies that designing of a general purpose optimizer is bootless. 
Independent of the fitness function, one cannot (without prior domain knowledge) 
successfully choose between two algorithms based on their success in a different set 
of problems. 


4 Evolutionary Computation 

Biological life appeared on the earth probably as a culmination of many chance 
events involving chemical and physical interactions of molecules. Ever since the 
appearance of unicellular life, progressively complex life-forms have gradually 
evolved by the process of genetic evolution. The process of genetic evolution has 
elements of adaptive learning built into it. Evolutionary Computing methodologies 
are mathematical models of natural evolution implemented on a computer and are 
the most versatile complex problem-solvers. These techniques operate on a set of in¬ 
dividuals or chromosomes (population) simultaneously. Each individual represents 
a potential solution to the problem being solved. The cardinality of population de¬ 
pends on the complexity of the problem. A chromosome is a sequence of genes 
(essentially the system parameters). Goodness of an individual is defined with a 
mathematical function, called fitness function (/: 5 -> R). Individuals having lower 
fitness value are slowly washed out by the dominant competitors. There is thus a nat¬ 
ural process of screening in the course of genetic evolution which finds expression 
in the Darwinian dictat of the survival of the fittest M- Chromosomal crossover 
and mutation produce new features in the chromosome. The evolution is a continu¬ 
ous process by which a species continuously strives to attain a genetic structure of 
the chromosomes that maximizes their probability of survival in a given environ¬ 
ment. The evolutionary algorithms continue their iterations to improve the fitness 
of the individuals until an fitness maximal solution is found. The way individuals 
in the population are distributed has a major influence on the search. In the pan- 
mictic model, there is no structure in the population, all individuals are potential 
partners. Each individual can interact with every other one in the population during 
the evolutionary process. However, it is possible to define some structure providing 
the algorithm with higher exploration capabilities of the search space with respect 
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to panmictic populations. There are two main canonical kinds of structured popu¬ 
lations, namely distributed or coarse-grained and cellular or fine-grained ElEl. 
The choice must be made carefully. 


Basic Ingredients of Evolutionary Computation 

In a diploid organism like the human being a pair of parental chromosomes form 
the chromosome-pair of the child, each chromosome containing millions of genes. 
One gene from the father and the corresponding gene from the mother constitute a 
gene-pair for the child, the gene-pair (genotype) determines the specihc attributes or 



Fig. 1 Flowchart for simple evolutionary computation 


characteristics called the phenotypes. Most gene values (alleles) are inherited by the 
child from the parents unaltered; but on rare occasions, one or more of them may 
undergo change, by a process called mutation. If such changes are benehcial for 
the individual for better adaptation to the environment, they have higher probabil¬ 
ities of survival, and have correspondingly higher probabilities for producing their 
offsprings. Over generations those benehcial changes tend to stay on while the non- 
benehcial ones tend to disappear. In Evolutionary Computation, new individuals are 
generated in the search space by applying certain genetic operators to the current 
population. The dominant individuals will mate more often creating descendants 
with similar or better htness in accordance with the natural selection process and 
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statistically moving toward more optimal places in the search space. A simple pic¬ 
ture of process steps in sequence for general Evolutionary Algorithms is shown in 
Eigure[T] The three basic components of EC are A) Selection, B) Crossover, C) Mu¬ 
tation. Each component can be realized in a number of different ways. Let us start 
by looking at different ways of implementing a selection process on an evolving 
population of potential solutions 


4.1 Selection 

Selection operator selects a proportion of the existing population with a specihed 
probability to create the basis for the next generation. The operator is designed to 
ensure that promising solutions of the population have a greater probability of be¬ 
ing selected for mating. The selection process is controlled by selection pressure, 
which is dehned as the degree to which the better individuals are favored. So, se¬ 
lection process is essentially biased towards the more ht individuals and drives the 
Evolutionary Computation (EC) to improve the population fitness over the succes¬ 
sive generations. Higher selection pressure can result in higher convergence rate, 
but with perhaps a higher chance of premature convergence. On the contrary, EC 
with low selection pressure will take unnecessarily longer time to find the optimal 
solution. In addition to the selection pressure, selection schemes should also try to 
preserve population diversity, to avoid premature convergence. The selection mech¬ 
anism consists of two steps, the selection probability calculation and the sampling 
algorithm. 


4.1.1 Generational replacement: 

Entire set of parents are replaced by their descendants or the n worst parents are 
replaced with n best offspring. 


4.1.2 Truncation Selection: 

A proportion (p) of the population are selected based on their fitness value and 
reproduced ^ times so that the population size is maintained. Less ht individuals 
are not given the opportunity to evolve. 


4.1.3 Roulette Wheel Sampling Algorithm: 

Individuals are selected according to their htness values Elm. The higher the 
htness, the higher the probability of being selected. Eor each individual C, in a pop¬ 
ulation P the selection probability {ps{Ci)) is given by 
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PsiCi) 


m) 


( 8 ) 


The population is then mapped onto a roulette wheel, where each chromosome C, is 
given a slot that proportionally corresponds to Ps(Ci). To select an individual, a ran¬ 
dom number is generated and the individual whose slot spans the random number is 
selected. The wheel is then spun N times, N being the cardinality of the population) 
to choose individuals until next generation mating pool is fully populated. 


4.1.4 Stochastic Universal Sampling: 

Population is represented by a pie chart, in which each zone is allocated for an in¬ 
dividual of the population. The areas of the zones are directly proportional to the 
fitness of the representative individuals, exactly as in roulette-wheel selection. But 
unlike roulette-wheel selection, here a single spin of the wheel is required to con¬ 
struct the mating pool. N (cardinality of the population) number of equally spaced 
pointers are put around the pie chart. The position of first pointer is generated ran¬ 
domly in the interval (0,1 /N). Individuals for mating pool are then selected by gen¬ 
erating the N pointers, starting with randomly generated first pointer, spaced by 
1/A^, and selecting the individuals whose fitness spans the positions of the pointers. 
Stochastic Universal Sampling provides low spreading over the desired distribution 
of individuals and is bias-free. This approach is better than roulette wheel, because 
it keeps the diversity and prevents best individual from dominating the population. 


4.1.5 Rank-based Selection: 

Individuals are ranked according to the fitness values in ascending order (least fit 
individual has rank=l, the fittest individual rank=A^). Rank based selection uses a 
function to map the indices of individuals in the sorted list to their selection proba¬ 
bilities. The performance of the selection scheme depends greatly on this mapping 
function. 

For linear rank-based selection, the bias is controlled by adjusting the selection 
pressure (SP), such that 2.0 >SP>\ .0, the expected sampling rate of the best indi¬ 
vidual being SP, while that of the worst individual is (2 — SP). The selection pressure 
for all other individuals in the population can be obtained by linear interpolation of 
the selection pressure according to the rank. The fitness value for individual is 
calculated as: 

F{i) = {2-SP)+2-(SP-\)-’^ (9) 

where F{i) and r, are the fitness and rank of the individual, respectively. Non¬ 
linear ranking permits higher selection pressures (N — 2 > SP > 1.0) than what is 
used in the linear ranking method. The corresponding fitness values can be com- 
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puted by using 



( 10 ) 


where X is the root of the polynomial equation 


{SP-N)-X^-^+SP-X^-^ + ... + SP-X + SP = 0 


( 11 ) 


Rank-based selection schemes help prevent premature convergence, but are compu¬ 
tationally expensive because of the need to sort populations, and slow convergence. 

4.1.6 Boltzmann Selection: 

A temperature like selection parameter controls this selection procedure. Initially 
the temperature has been kept high and then gradually lowered (selection pressure 
increases gradually) throughout the EA run. The probability of accepting an indi¬ 
vidual k: 



( 12 ) 


where and Fa are the fitness of kf^ individual and average fitness of the population 
respectively. Kb is Boltzmann constant and T is the temperature. 

4.1.7 Tournament Selection: 

A tournament is held among n randomly picked competitors from the population. 
The individual with the best fitness value of those random n tournament competitors 
(winner of the tournament) is then copied into the mating pool. A tournament size 
of n = 2 is commonly used in practice. Selection pressure can be controlled by con¬ 
trolling the tournament size n. Tournament selection is superior to the proportionate 
selection as the latter is found to be significantly slower in terms of convergence 
time. Linear ranking and stochastic binary tournament selection have identical per¬ 
formance in expectation, but binary tournament selection is preferred because of its 
more efficient time complexity 1201 . 

4.1.8 Sexual Selection: 

Crossover takes place only between chromosomes of opposite sex. The sex of the 
individuals in the current generation are determined either randomly or based on 
some specific criterion of the individuals. Generally all the females will get to re- 
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produce only once regardless of their fitness level to facilitate exploration of the 
search space. The male selection is fitness biased. In feminine selection, the female 
chooses one of the males to mate based on a problem dependent attraction function 

MM. 


4.1.9 Clonal Selection: 

A proportion of the fitter individuals are selected for cloning. Each of these indi¬ 
viduals receive a number of copies proportional to its position in the ranking. The 
clones then undergo the maturation process. A given individual and its maturated 
clones forms a subpopulation and the best of each subpopulation is allowed to pass 
to the next generation 1271 . 


4.2 Crossover 


It is designed for sharing information between individuals by swapping or intermin¬ 
gling the genetic materials of two randomly chosen parent chromosomes, with the 
possibility that good chromosomes may generate promising descendants. A user de- 


Parent strings 



0 B 

^[d 

E 

F 

G [h 

I 


[a 

B 
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® @(D(D 


Crossover site (randomly chosen) 


Fig. 2 Simple Crossover involves exchanging genetic information, beyond a randomly chosen 
cross site, by swapping the gene values within the parent’s chromosome. Traditionally chromo¬ 
somes are crossed at a single point. However, some problems could benefit from using multiple 
crossover points. 


fined crossover probability (pc) determines the crossover frequency, i.e., how often 
crossover will take place. A general crossover scheme is shown in Eigurej^ In the 
discussion of different crossover techniques that follows, we have represented the 
parent and offspring, of cardinality d, in a population as; 


Parent 

Pj{PvPi,-,Pd) 


Offspring 
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4.2.1 Simple crossover; 


A single crossover point / S {2, 3, — 1} is randomly chosen (say, i), where 

the two mating chromosomes are cut and the sections after the cuts are swapped 
to create two new individuals. N point crossover is a generalization of this simple 
crossover applied to N different segments. 


<^1(01, 

Oiioi, 


02,-,0d) = {p{,Pi,Pv ■■■,pi,\ phi, ■■■,Pd-vPd) 
02,-,0d) = {p\,P2,Ph ph'^ Pi+V ■■■,Pd-vPd) 


4.2.2 Cut and splice crossover: 

This operator is an extension of simple crossover and designed specifically for clus¬ 
ter geometry optimization problems (see section Metal Clusters & Nanoalloys). 


4.2.3 Uniform crossover: 

It is assumed to be more disruptive in nature than simple crossover. So a lower 
crossover rate (say, 0.50) is used. Descendants are created by swapping the genes of 
two parent chromosomes with a predefined probability M- 


Ipj, if n>0.5 
I fh otherwise 


n e [0,1] 


uniform random number 


02(01,oi...,o^) = ( p\, 


(13) 
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4.2.4 Cell crossover: 

Cells (a set of genes) of two random parent individuals are swapped with predefined 
probability to produce two descendants ll25l . 
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4.2.5 Heuristic crossover; 


A single offspring is generated by linear extrapolation of two parent individual using 

It is 


equation 14 that may lie outside of the range of the two parent vectors 
also possible for this operator to generate an unfeasible offspring vector. In such a 
case another offspring created with different choice of r, r being a uniform random 
variable belonging to the interval [0,1]. 


{oi} = {pi +ripi-pi)} subject to F{Pj) > F{Pk) 


(14) 


4.2.6 Simplex crossover; 

Simplex is constructed with randomly selected N parents from the population. Cen¬ 
troid (C) of the simplex is calculated by 


N p 
!=0 

The simplex is then expanded by a small degree e 

Oj=C + e-{Pj-C), je {1,2,-, N} (16) 

and within this expanded simplex, one or more new individuals are sampled ll27l . 



I Single point Arithmetic Crossover 

randomly chosen crossover point (s) 

/ I i i 



Multi point Arithmetic Crossover 


Fig. 3 Arithmatic Crossover involves intermingling genetic information, beyond a randomly cho¬ 
sen cross site and creates descendants that are the weighted arithmetic mean of two parents. For 
simplicity all the gene values in the parent strings are shown by same color. 
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4.2.7 Linear crossover; 

Three offspring are generated (one in the exploitation zone and two in the explo¬ 
ration zone), with the two most promising solutions being chosen as descendants by 
means of fitness evaluation EEl . 



(17) 


4.2.8 Arithmetic crossover; 

Like simple crossover, chromosomes are crossed at a single or multiple crossover 
points. All positions arithmetic crossover are also useful. It was introduced to pre¬ 
vent invalid crossover and to increase searching space and convergence speed. Two 
complimentary linear combinations of the parents are generated on the basis of the 
arithmetic mean . 


Oi = aPk + {l-a)Pj 

02 = aPj + ll-a)Pk 


(18) 


Arithmetic crossover operates in a way that every gene in descendants is a con¬ 
vex combination of the cooresponding genes from the two parents, a is constant 
in uniform arithmetical crossover or it may vary with generations in non-uniform 
arithmetical crossover. Figure [^schematic ally shows arithmatic crossover. 

4.2.9 BLX-a crossover; 

To expand the range of arithmetic crossover Eshelman and Schaffer suggested blend 
crossover ll29l . which generates a single offspring by blending two floating point 
parent vectors as follows: 


oj = 91 {{LBi -a-1), {UBi + a-I)) 
where, UBi = max {pj , pf), LBi = min (p/, pf), I = UBi — LBi 


(19) 


3i{a,b) is a function generating a uniform random number in the range {a,b). The 
user-defined parameter a is usually set to a value of 0.5. Setting a = 0.0 and 0.25 
will give flat crossover OOl and extended intermediate crossover ED, respectively. 

4.2.10 Similarity crossover; 

In arithmetic crossover, there has always been a problem in selecting a perfect 
weighting factor. The weight factor (a) in arithmetic crossover (equation [TSl l is 
replaced by similarity measurement (5 = between parents 
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EN=~{ {p\ -pT? + {p ‘2 - P'^f + --- + {p'd-p7?f (20) 

In addition to these, there are other variants of crossover schemes such as parent¬ 
centric BLX-a crossover ll3^ . Fuzzy recombination ll34ll . SBX lIlSll . PCX IMIl . 
XLM IIJTII . Laplace crossover ll38l . differential evolution ll39l . partition ll40l l4Tll . 
linear BGA 1421, UNDX ll43l . fuzzy connectives based ll44l . direction-based ll45ll . 
multiple crossover P6l , Ring crossover ill, hybrid crossovers 1481 are available in 
the literature. 


4.3 Mutation 

The mutation operator perturbs one or more components (genes) of a selected chro¬ 
mosome, regulated by a predefined mutation probability, so as to increase en¬ 
tropy or to decrease the mutual information in the population. Mutation simply re¬ 
stores lost information or import unexplored genetic material into the population in 
order to distribute solutions widely across the search space and thus avoid premature 
convergence. 

Before mutation After mutation 


(D©©^E 


Mutation site (randomly chosen) 

Fig. 4 Mutation causes small alterations at randomly chosen gene(s) (with a predefined mutation 
probability) in an individual within available range of possible values. 


The performance of any evolutionary search algorithm heavily depends on the 
choice of its main control parameters: population size, mutation rate, and recombi¬ 
nation rate and among these mutation rate is the most sensitive control parameter. 
Mutation-only-EA is possible, but, crossover-only-EA would not work with finite 
population size. There should be always a proper balance between exploration (dis¬ 
covering promising areas in the search space) and exploitation (optimizing within a 
promising area already found) abilities of a search algorithm. Crossover exploits ex¬ 
isting genetic material, while mutation explores for newer material. Mutation does 
its best to avoid premature convergence and explore newer areas of search space. 
High mutation rate increases the probability of searching over large areas of search 
space; however, it prevents the population from converging to an optimum solution. 
On the other hand, too small a mutation rate may result in premature convergence. 
The best value of mutation rate is strongly problem specific and depends on the na¬ 
ture and implementation of the algorithm. There is no universal best mutation rate 
for most of the real world problems. The crossover and mutation rates (probabil¬ 
ity) are critical to the success of any EA, and are usually determined by means of 
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trial-and-eiTor. Automating the parameter settings i.e., making the search algorithm 
self-adaptive has been one of the most fashionable area of current interest among 
the EA community. A general mutation scheme is shown pictorially in Eigure 
We shall use specihc notations for an individual before and after mutation, right 
through the discussions of different mutation schemes. The chromosome before and 
after mutation are expressed as follows; 

Before mutation After mutation 

X = {xi,x2,--- ,xi,--- ,xd-\,xd) Y = {yi,y2,--- ,yi,-" ,yd-\,yd) 


4.3.1 Uniform mutation: 

Genes in an individual are selected with a predehned probability (pm) and replaced 
with a uniform random value according to ll28ll . 

^ 1 3i{LBi,UBi), ifp,„ > r,-, r,- S (0,1) 

]xi, otherwise 

3i{a,b) is a function generating a uniform random number in the range {a,b), LBi 
and UBi represents lower bound and upper bound to the gene. Uniform mutation 
tends to forestall premature convergence. In the neighborhood of the optimal region, 
however, uniform mutation rate may disrupt the hne-tuning of the search to the 
optimal point by inducing random oscillations. 


4.3.2 Nonuniform mutation: 

Due to its dynamical nature, non-uniform mutation 12^ is one of the most com¬ 
monly used mutation scheme in real coded EC. Gene(s), selected with a predehned 
mutation probability (pm), are mutated as 


Xi + A{t,UBi-Xi), ifr > 0.5 
Xi — A (f, Xi — LBi), otherwise 

A(f,z) =z^l —^ 

where r G [0,1] is a random number, and t is the generation number. The user dehned 
parameter b takes care of the non-uniform nature of the mutation step sizes. 



( 22 ) 

(23) 


4.3.3 Gaussian mutation: 

A noise with a mean of 0 and a standard deviation of 0.1 times the maximum value 
of the gene is added to the gene-value. A new individual is obtained by adding a 
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random value to each element of the selected parent. 


_ ( Xi+N{0,ai) Gaussian mutation ESl 

|x, -|-C(0,T,) Cauchy mutation ll4^ 

The mutation function that describes the transformation of an element x into y l49l 


G= J Gaussian (25) 

“y 

T is the scale parameter. Ewfeevy-type mutation 1501 im, Real Number Creep 1521 
are some variants of Gaussian mutation. The scheme of Real Number Creep; 


yi = 


Xi+N{Q,af) \fi = k 
Xi otherwise 


UBi-LBj 
= - 1000 - 


k = Uniform[l,(f] 


(26) 

(27) 


Based on the PBX-a recombination operator and Gaussian distribution function, 
Dorronsoro et al. designed GPBX— a l53l mutation operator. 


4.3.4 Dynamic Random Mutation: 

The scheme enables the mutation size to be adjusted dynamically l54ll . 

yi=Xi+s-4>o-{UB-LB) (28) 

where is a random perturbation vector within the interval of (0,1]. The mutation 
step size is dynamically tuned by 

(29) 

\ ^max J 

The parameter b > 0 controls the decay rate of s and t^ax denotes the maximum 
number of generations over which the search takes place. The parameter b gov¬ 
erns the shape of the allowable mutation region. The mutation range dynamically 
decreases as the number of generation increases. The idea is adopted from the an¬ 
nealing procedure in metallurgy. 


4.3.5 Directed Random Mutation: 


The value of a randomly selected gene is modified according to 


J Xi ±Am- r, if Pm < rand (0,1) 
1^ Xi, otherwise 


(30) 
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where r is a random number in the interval 0 < r < 1. The + or — sign is chosen 
with a probability of 0.5. A„, is the mutation intensity, dynamically adjusted on the 
basis of the degree of acceptability of mutation over a number of past generations 

isiisa. 


( {I + Kg), if Naccept < Nt ower) 

if Naccept ^ Supper) 


where Pg is a Gaussian random number in the range (0,1) and Naccept is the num¬ 
ber of accepted generation over past N generations. Ntower and Nupper are the user 
defined lower and upper bound to Naccept- 


4.3.6 Wavelet mutation: 

A continuous-time function is called a “mother wavelet” if it satisfies the 

following properties; 


^+oo 

I \j/{co)dco ^0 

(32) 

.+O0 

1 \\j/{a))\ d(0 < oo 

(33) 


One example of such mother wavelet is Morlet wavelet: 

_(£?■ 

\l/{co) = e ^cos{5co) (34) 

In order to control the magnitude and the position of yf{co) , V4.t)(®) is defined as 

Va,bi^) = (35) 

y'a \ a J 

where a is the dilation parameter and b is the translation parameter. Ling et al. 
Il57l . based on the wavelet theory, proposed wavelet mutation where each gene in 
a chromosome have a chance to mutate according to a preset mutation probabilty 
Pm € [0,1]. The mutated gene is represented by 

^ l'xi + i/Pa,o((0)(l7B-Xi) ifV4,o(®)>0 

lxt + l//a,o((0) (xt-LB) if Vet,c((0) < 0 ^ 

For Morlet wavelet o G [—2.5,2.5] is randomly generated. 


4.3.7 BGA mutation: 

The mutated gene y, is computed according to ED 
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yi = Xi ± ranget ■ 5 (37) 

The + or — sign is chosen with a probability of 0.5, ranget is normally set to 0.1 • 
{UBi — LBi) and 5 is computed as follows. 

d 

5 = £ ai2-‘ Ai e 0,1 (38) 

i=i 

a, is randomly generated and takes the value 1 with a probability Ps = d- The 
scheme ensures that at least one variable will be mutated in an individual. Discrete 
and continuous modal mutations Il58l are generalization of EGA mutation. They 
only differ in how 5 is computed. 


4.4 Supplement to Selection: Elitism 

Sometimes crossover and mutation operations destroy potentially valuable genetic 
information which has evolved during the search. There is no guarantee whether EA 
will re-discover these lost information or not. If EA rediscovers those previously 
discarded information, it does so at the cost of some unnecessary time and compu¬ 
tational power. To prevent this, the concept of elitism is introduced. Elitism ensure 
that a small proportion (typically 1 or 2) of the fittest individuals are passed onto 
the new generation so that the valuable information remains intact within the popu¬ 
lation and thus reduces the genetic drift. Higher proportion (causes rise in selection 
pressure) may lead to premature convergence. So the degree of elitism should be 
adjusted carefully. Elitism can speed up the performance of the GA significantly. 

Apart from these operators, initial population also has major influence on the 
convergence of the algorithms. It is difficult to hit the global optimum in large sys¬ 
tems starting from a fully randomized initial population as it leads to nearly iden¬ 
tical glassy structures that have similar (high) energies and low degree of order 
1591 . Often random but symmetric structures may be helpful as initial guess for sys¬ 
tems with large dimensions llbOll . Presence of constraints often make an optimization 
problem difficult to solve. Several techniques such as adding penalty term to the ob¬ 
jective function (the penalty increases with constraint violation), adjustment of the 
crossover or mutation operation, inducing instabilities to the infeasible solutions of 
the populations without modifying the objective function, have been proposed llbTII . 


4.5 Evolutionary Computing Methodologies 

Let us consider several variants of EC methodologies that have been in use. 


20 



K Sarkar 


Pure and Hybrid Evolutionary Computing in Global Optimization 


4.5.1 Random Mutation Hill Climbing: 

It is the simplest algorithm in the class of Evolutionary Computing methodologies. 
The algorithm works with one randomly generated string containing all the system 
parameters, (say, Cartesian position variables for geometry optimization) required 
for optimization and it involves mutation, as the only evolutionary process in order 
to climbing up the htness landscape. Thus the underlying principle is to randomly 
choose a solution in the neighborhood of the current solution by means of mutating 
the string and retain this new solution if it improves the fitness function. A modihed 
Random Mutation Hill Climbing was proposed and implemented by Sarkar et al. 
Il56l|62l|63]l6l. The modification lies in additional built-in features for enforcing 
adaptive control on all the parameters of the search heuristic. 


4.5.2 Genetic Algorithms: 

GAs are a class of stochastic heuristic search methods based on simplifications of 
evolutionary processes observed in Nature. They maintain a finite population size 
and are typically good at both the exploration and exploitation of the search space. 
They are invoked particularly in large, complex, and poorly understood problem 
domain. The search begins with randomly initializing a population of individuals 
followed by selection operation to create a mating pool. The mating process, imple¬ 
mented by crossing over the genetic material from two parents to form the genetic 
material for two new solutions. Random mutation(s) is(are) applied to promote di¬ 
versity. If the new solutions are better than those in the population, the individuals 
in the population are replaced by the new solutions. The key idea is to maintain a 
population of chromosomes, that evolves over time through a process of competi¬ 
tion and controlled variation. GA have had a great measure of success in search and 
optimization problems due to their ability to exploit the information accumulated 
about an initially unknown search space in order to bias subsequent searches into 
useful subspaces. The population scalability of performance in evolutionary com¬ 
putation has been an issue pursued since 1990s. The population sizing equation of 
Goldberg ll65l . Miihlenbein’s study on minimal population size that can ensure con¬ 
vergence to the minimum Il66l . the work relating to adaptive control of population 
size by Arabas et al ill are the early examples of interest in this area. However, 
GA is trivially parallelized and its parallel efficiency could often outweigh the higher 
computational demand in every generation . 


4.5.3 Evolutionary Programming: 

Both Evolutionary Programming (EP) and Evolution Strategy (ES) are known as 
phenotypic algorithms, whereas the GA are genotypic algorithms. Phenotypic al¬ 
gorithms operate directly on the parameters of the system itself, whereas genotypic 
algorithms operate on chromosomes or individuals representing the system. In re- 
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cent work however the distinction has become blurred as ideas from the ES and EP 
have been applied also to strings or chromosomes that constitutes a population in 
GA. The primary difference between EP and the other approaches (GA, GP, and 
ES) is that no exchange of material between individuals in the population is made. 
Evolution is caused solely by appropriate mutation operators. 

4.5.4 Evolution Strategy; 

The earliest form of ES is based on populations that consist of two competing in¬ 
dividuals and the evolution process consists of the two operations mutation and 
selection. Only the fittest of the two individuals is allowed to produce descendants 
in the following generation. Later, general versions of the algorithm were devel¬ 
oped through increasing the number of members in a population (/r-parents and 
A-descendants), among which the most common strategies are ()i,A) and (jj. -f A). 
In (/X, A) the parents are simply replaced with the descendants and in (jj. -f A), the 
parents compete with their children. 

4.5.5 Differential Evolution; 

DE is a population based technique 1^17011711 . The members of the population are 
chosen randomly to start with (Xi,X 2 , ■ ■ ■ ,X„). In a Differential Evolution (DE), mu¬ 
tation plays the dominant role. Mutation is used to produce new individual by adding 
the vector-difference between two randomly selected members (say, Aj. and A/) of 
the population to another (say, A„,) producing a trial vector V,, where. 


Vi=X„,+f{Xk-Xi) 


(39) 


/ being the mutation intensity. A total of n such individuals are generated (Vi, V 2 , 
• • •, V„). A crossover operator is then used to produce new candidate t/,- from the 
trial vector V, by setting the component of V) as follows: 


— E/y if r(y) ^ Pc or j — j random 
= A,y, Otherwise 


(40) 

(41) 


The candidate is evaluated, and the better of the current individual and the candidate 
solution is selected. 

4.5.6 Genetic Programming; 

GP is an extension of the GA in which the structures in the population are not fixed- 
length strings that encode candidate solutions to a problem, but programs expressed 
as syntax trees. So essentially Genetic Programming is a method to evolve computer 
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programs. The genetic operators allow the syntax of resulting expressions to be 
preserved. 


4.5.7 Memetic Algorithms: 

MAs are extensions of EA that repeatedly apply a local-search algorithm to exploit 
the search space in between generations to reduce the likelihood of the premature 
convergence. 


4.5.8 Particle Swarm Optimization: 

PSO was introduced by Kennedy and Eberhart fflX as a gradient free evolution¬ 
ary computation technique that utilizes the ‘collective intelligence’ of a coopera¬ 
tive swarm of particles. Each particle representing a candidate solution is randomly 
assigned variables to be optimized and then allowed to “fly” with random veloc¬ 
ities < V < Vfnax) in the search space. The task is to reach a point in the N- 
dimensional search space where the htness has the maximum value and which is 
hopefully the optimal solution as well. Each particle keeps track of its coordinates 
in the search space, associated with the previous best fitness value ipbest)- The over¬ 
all best htness value that the swarm has achieved so far is stored as gbest- At f -f 1 
generation, velocity (v) and position (x) of particle are updated based on its cur¬ 
rent velocity (v'f), the distance from its previous best position (pj), and the distance 
from the global best position (g'j) according to 

= a)v\+Cl(t>l{p'i-^^i}+C2(t>2{g'i-^i} (42) 

x'+*=4-fv'+' (43) 

The updating rule applies to all the components of velocities and positions. Ci,C 2 
are the acceleration coefficients usually set to a value of 2.05, and are uni¬ 
formly distributed random numbers within the interval of (0,1). The inertia weight 
(0 decreases linearly from (O — 0.9 to 0) = 0.4 over the whole run to control the mo¬ 
mentum of the particles. The coefficient Ci is called particle’s self-conhdence due 
to its a contribution towards the self-exploration of a particle, whereas C 2 is called 
swarm conhdence due to its contribution towards motion of the particles in global 
direction, which takes into the collective wisdom GS. The velocities are assigned 
cutoffs Vmin and Vmax to prevent the swarm from disintegrating. The coordinates 
evolve continuously, and have no cutoffs. 
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5 Monte Carlo Algorithms 

The main principles of Monte Carlo algorithm, named after a famous casino city in 
Monaco, are ergodicity and repeated random sampling to find the optimal solution. 
It was originally practiced under more generic names such as statistical sampling. 
In the 1940s, physicists working on nuclear weapons projects at the Los Alamos 
National Laboratory coined the present name. The random search, direct Monte 
Carlo sampling, the random hill climbing methods. Simulated Annealing, Quantum 
Annealing and Parallel Tempering are some examples that belongs to the class of 
Monte Carlo algorithms. 

Simulated Annealing; 

The term “annealing” is commonly used in metallurgy to define a process in which 
the metal is heated to a high temperature, followed by slow and controlled cooling. 
If molten metal is quickly cooled, it does not reach the perfect crystalline state of 
minimum free energy. The essence of the process lies in slow cooling which allows, 
time for redistribution of the atoms as they slowly lose mobility. Metropolis et al. 
Il74l proposed an algorithm for finding the equilibrium configuration of a collection 
of interacting atoms at a given temperature. The transition to a new state (j) from a 
state (i) is accepted with the probability 



(44) 


otherwise 


Kirkpatrick et al. US introduced the concept of Simulated annealing (SA) in com¬ 
binatorial optimization following Metrpolis’s idea. In SA, energy is regarded as the 
cost function and a set of parameters {x, } is used to define the configuration at some 
effective temperature. The target is to find the global minimum of the cost function 
using temperature as the control parameter. The SA first allows the system to ‘melt’ 
at a high temperature and the temperature is lowered in small steps allowing it to 
spend sufficient time at each temperature until the system freezes. At each step of 
temperature, the sampling must be carried out for a sufficiently long time for the 
system to reach a steady state. The sequence of temperatures and the number of 
samplings allowed to reach ‘equilibrium’ at each temperature constitute an “anneal¬ 
ing schedule”. At each reconfiguring step perturbations are applied to generate a 
new point in the search space. The new point is accepted based on the Metropolis 
algorithm mentioned above (equation |4^. Thus the algorithm differs from the Hill- 
Climbing method in the criterion used to replace the existing solution with a new 
one. The advantage of S A is the built-in local optimization scheme, since in the limit 
of the vanishing control parameter (T), the algorithm is a standard downhill simplex 
method. A performance comparison with LAs shows different strengths and weak¬ 
nesses. While the short-range exploitation of SA based algorithms is generally con¬ 
sidered to be better than that of pure LAs, LAs generally provide a better long-range 
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exploration. Employing hybrid EAs with enabled local optimization steps in general 
makes the EAs at least on par with SA-based methods. 

An inverse approach (Reverse Monte Carlo) to the classical Monte Carlo algo¬ 
rithms successfully reproduces many salient features of glass and metallic glass 
systems obtained from molecular dynamics simulation 17611751. Along with the 
success in different optimization problems SA has also been effectively used in 
quantum mechanical calculations, especially in the optimization of parameters of 
wave functions, constructing reaction paths for transformations on model potential 
energy surfaces and for rearrangement in Lennard Jones clusters iTTTllTSll . 

We have so far summarized the important aspects of evolutionary computing along 
with some relevant natural computing methods. We hope the readers will be able to 
exploit the information to form a well founded basis to choose appropriate schemes 
which maintain a good balance between exploration and exploitation for specific 
problems. In the next section we will provide an overview of applications of evo¬ 
lutionary computing and its hybrid methods in the general context of computing 
electronic structure of atoms, molecules and clusters. 


6 Evolutionary Computing in Electronic Structure Calculation 

The applications of EC to problems of calculation of electronic structure of atoms, 
molecules, clusters and crystals have been interesting and numerous. We will review 
the applications under these categories: 

• Solution of Schrodinger equation (SE) by EC. 

• Locating minima on the complex potential energy surfaces of molecules, clusters, 
crystal structure prediction. 


6.1 GA and Electronic Structure Calculation 

The first attempt to marshal the power of GA to tackle the problem of electronic 
structure calculation can be traced to a paper by Zeiri et al. 1791 on the application 
of GA to the calculation of bound state within the framework of a local density 
approximation almost 30 years after GA came into being. GA was invoked for di¬ 
rectly solving the time-independent SE ISOl by Chaudhury and Bhattacharyya in 
1998. They tested the workability of their algorithm in finding the ground state of a 
particle moving in screened coulomb potential and that of an oscillator with quartic 
anharmonicity. They went on to examine the applicability of their method to the 
problem of solving the inhomogeneous differential equations of the RSPT with spe¬ 
cial reference to the ground states of two-electron atoms. Nakanishi and Sugawara 
(2000) proposed a ANN based way of solving the SE in which the wavefunction was 
represented on a perceptron type NN, and the weights and the biases were optimized 
by a micro-GA ISTl . 
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Saha and Bhattacharyya (2001) proposed a stable and generalizable basis set 
free strategy of solving the SE by GA ||8^ . Test calculations on the ground and 
excited states of H-atom, coupled harmonic oscillators, and symmetric double well, 
demonstrated the robustness and workability of the algorithm proposed. Saha et al. 
(2003) parallelized ll8^ the fitness evaluation steps and the parallel GA method was 
successfully invoked to solve an exactly solvable coupled oscillator problem and 
the H-atom problem. The kinetic energy was computed by numerical EET while 
the potential energy was obtained by quadrature. Saha and Bhattacharyya (2004) 
further explored GA in solving the radial SE for the ground state of two electron 
atoms and ions ll84l . The approach is basis set free and the results were better than 
HE results. The same authors extended the theory to calculate the ground state of 
Hydrogen molecule ion {H^)- In one realization, the internuclear distance (R) was 
allowed to evolve along with amplitudes of the electronic wave function at different 
points on a 3-d uniform grid. The equilibrium internuclear separation and the energy 
was obtained in a single run. Sahin et al. (2006) explored the GA based approach 
for computing the energy levels of hydrogenic impurity in a quantum dot and then 
went on to further extend the method for the self-consistent electronic structure 
calculation of many electron quantum dots |[85l . 

GA has been deftly used to compute the molecular electronic structure of doped 
as well as undoped polythiophene and polyselenophene oligomers ll5^ [86l within 
the framework of a modified SSH method. The density based approach in which the 
GA operators act only on the molecular geometry strings (arrays of the geometrical 
parameters) that uniquely define the corresponding fixed nuclei Hamiltonian H{R). 
H{R) in turn, defines a unitary transformation U{R) which transforms a trial density 
P{Ro) into a new density P(/?) = U{R)^P{Ro)U(R) at the new genetically modified 
geometry. Thus along with the geometry strings, the density P{R) is forced to co¬ 
evolve, till [H{R),P{R)] = 0. The technique appears to be sufficiently flexible and 
could find use in ab-initio calculations of molecular electronic structure. 


6.2 GA in the Exploration of PES of Clusters 

Traditional optimization methods perform well on strongly convex response sur¬ 
faces, whereas EAs are particularly appropriate for most of the real-world NP-hard 
problems in which the objective is multimodal and expensive to evaluate, the gra¬ 
dient of the objective function is not analytically definable, the problem domain is 
large in dimension and/or multiple linear and non-linear constraints have to be ap¬ 
plied simultaneously. Traditional hard computing methods are occasionally fused 
with the evolutionary methodologies to develop computationally intelligent hybrid 
systems with moderate computational cost. 

Atomic or molecular clusters are aggregates of similar or dissimilar particles and 
generally have physical and chemical properties different from the bulk. Geometry 
optimization of atomic or molecular clusters belongs to the class of nondetermin- 
istic polynomial hard (NP-hard) problems. Exhaustive search throughout the PES 
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in any reasonable amount of time is almost an impossible task. By the early 1990s 
EAs were applied to these classes of problems. Deaven et al. lISTlI applied GA to 
find fullerene structures up to Cgo starting from random atomic coordinates. Empir¬ 
ical Brenner potential was used to model the PES. Dugan et al. If88l had proposed 
Monte Carlo type local optimization between GA steps for moderately sized carbon 
clusters. Chaudhury et al. Il89ll have used EA methods on empirical potentials to ob¬ 
tain optimized structures of halide ions in water clusters, which they then coupled 
with quantum-chemical (AMI) calculations for simulation of vibrational spectra. 
By suitably defining the objective functions, EC methods can also be used to locate 
additional features on the PES rather than just low-energy local and global minima. 
Chaudhury et al. Il90l 1^ have implemented methods for finding first-order saddle 
points and reaction paths, on PES of Lennard-Jones (LJ) clusters up to n=30. 

Clusters may have properties which are different from those of discrete molecules 
or bulk matter: for example, some metals (e.g. palladium) which are non-magnetic 
in the solid state, are magnetic, with relatively high local magnetic moments in 
both neutral and anionic PdAj clusters ll92l . The Coulomb explosion of large rare- 
gas clusters release huge energy (approaches the energy of nuclear processes). At 
very low temperatures (< 2K for ^He), He clusters may form a superfluid droplet. 
Metal nanoparticles play a crucial role as heterogeneous catalysts in petrochemi¬ 
cal, pharmaceutical and clean energy sectors. Novel and interesting nanostructures, 
having useful chemical and physical properties could be obtained by controlling the 
nanoparticle size, composition, surface site preferences, and degree of segregation 
between the metal constituents ll^ . Therefore, computing ground state configu¬ 
ration of atomic and molecular clusters is a challenging task from the theoretical 
perspective. The number of local minima of the clusters rises exponentially with the 
growth in cluster size. It increases further with composition heterogeneity because 
structures more complex than their homogeneous counterparts are possible. It is due 
to the existence of isomers with the same geometry and composition but showing 
different distributions of the constituent particles that the optimization task becomes 
notoriously difficult ||95l [TOOl [M] |96l |92l IHHE |94l [Ml- For a cluster made of 
Ha and ng atoms, with nA + rig — N, one finds that a single geometrical isomer can 
have , different homotopes. Various attempts are there in the literature to And 
the global minimum energy structure of these systems, as for example the dynamic 
lattice searching method 11031 . basin-hopping (BH) approach lll04lll05j|106l[T07ll . 
adaptive immune optimization algorithm (AIOA) 0108111091II10111111 and evolu¬ 
tionary algorithms ill HIE El EE El EE ES EE EE El EE EE 

1 12411125ll 12611127111281 . Hybrid approaches combining evolutionary algorithm and 
a local search procedure, have been increasingly used to handle these problems. We 
review recent work in the area in what follows. 


6.2.1 LJ and Morse Clusters 

Ab initio quantum mechanical calculation for finding the global minimum geometry 
of a cluster of atoms or molecules is often extremely expensive II129II1301 . Compu- 
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tational cost increases with the accuracy of the calculation and size of the system. 
Generally empirical or semiempirical pair potentials are employed for large systems. 
LJ and Morse functions are the two most widely used and relatively simple pairwise 
additive potential models describing interactions between atoms which have been 
regularly used as benchmark to assess the performance of global search methods for 
cluster geometry optimization. The LJ and Morse potentials are respectively given 
by 



(45) 



(46) 


1=1 i>i 


where stands for the Euclidean distance between atoms i and j, rg is the equi¬ 
librium bond length, e is the well depth of the pair-potential, (7 is the separation 
at which the pair-potential between the atoms goes through zero and p is the in¬ 
teraction range scaling parameter. Eor the prediction of the ground state structure 
of crystals, the USPEX method ifTIll fT^ fT^ fml fml fT4Ql fml 

has turned out to be extremely powerful and already guided materials scientists in 
finding interesting and unexpected crystal structures. Zhu et al. 11411 enhanced the 
efficiency of USPEX method by designing additional variation operators and con¬ 
straints for partially or completely fixed molecules. Goedecker et al. Ill42lll43lll44l 
offered improved minima hopping with a softening method and a stronger feedback 
mechanism to predict structures of homoatomic and binary clusters with LJ inter¬ 
action as well as structures of silicon and gold clusters described by force fields. 
Cheng et al. Il97ll employed funnel hopping algorithm coupled with GA to locate 
the putative global minima of the LJ clusters and the Morse clusters up to N = 160. 
Wales et al. showed that the use of approximate symmetry provides a more produc¬ 
tive way to explore the configuration space and substantially improves the efficacy 
of global optimization for the atomic clusters 014511146111471 Eroltsov et al. 
il49l developed the “cut and splice” GA augmented with twinning mutation moves 
for the structural optimization of both a series of single-funnel LJ clusters up to 
70 atoms and the double-funnel LJ 38 cluster. Binary Lennard Jones (BLJ) cluster 
optimization is an even more challenging problem from the point of view of com¬ 
binatorial complexity. They are also interesting because catalytic properties of such 
clusters depend on the composition and structure. Comparatively little work has 
been done on the mixed clusters II1501115111152II as the search space is dramati¬ 
cally enlarged with the inclusion of more atom types. Mixed LJ clusters of widely 
varying compositions offer highly interesting additional perspectives on how vari¬ 
ations in preferred structures emerge. Hybrid approaches like combining a global 
optimizer with a local search procedure have proved to be useful for mixed clusters 
015311154lll55lfT5^ll571ll58lll59l . 17 new putative global minima for BLJ clusters 
in the size range of 90-100 particles have been predicted by coupling hidden-force 
algorithm and non-Markovian parallel Monte Carlo search 11551 . Dor et al., pro- 
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posed a multi-swarm based algorithm called PSO-2S 115611 and tested it success¬ 
fully on a number of benchmark functions. It uses charged particles in a partitioned 
search space and applies a repulsion heuristic on particles to increase the diver¬ 
sity for solving both unimodal and multimodal problems. Using GA-based global 
structure optimization framework and new set of fitted parameters for LJ potential 
from high-end ab initio calculations Hartke et al. optimized strongly mixed binary 
to quinary rare gas clusters 0157111581 . An average offspring method 01421 designed 
for cluster structure prediction was shown to perform better in systems with compact 
optimal structures. Two individuals are randomly selected and for each atom of the 
first cluster, the closest lying atom of the other cluster is identified. The correspond¬ 
ing atom of the child is placed randomly on the connecting line between the two 
parent atoms. The randomness of this operator is necessary to prevent the algorithm 
from producing a lot of identical offsprings. Leitao et al. II1601 claimed that they ap¬ 
plied for the first time an Island Model to the optimization of Morse clusters ranging 
from 41 to 80 atoms, combined with a hybrid steady-state evolutionary algorithm 
and a local optimization method. The performance was slightly more robust than 
that of the sequential approach. Dieterich et al. II15811161111621 designed EA pro¬ 
gram suite OGOLEM for structure optimization of mixed clusters. The OGOLEM 
framework provides both an MPI fronted for MPP parallelization and a threading 
fronted for SMP parallelization omitting unnecessary MPI overhead. They demon¬ 
strated the possibility to design molecules with targetted properties in the area of 
photochemistry using the OGOLEM framework III 6 II . There are prescriptions for 
the inclusion of Taboo-search features 0163111641 into the evolutionary algorithms 
which might help in reducing the amount of time spent in local optimizations re¬ 
discovering already known minima. Daskin et al. 11651 presented Group Leaders 
optimization algorithm in which the influence of the leaders in social groups is the 
inspiration for the evolutionary technique. The method is applied to locate the geo¬ 
metric structures of LJ clusters as well as to the quantum circuit design problems. 
Deep et al. II1661 made an attempt to solve LJ problem by incorporating a multi-orbit 
dynamic neighborhood topology in PSO. In multi-orbit topology, the swarm has het¬ 
erogeneous connectivity with some subsets of the swarm strongly connected while 
with the others are relatively isolated. This heterogeneity of connections balances 
the exploration-exploitation trade-off in the swarm. The dynamic neighborhoods 
topology helps avoid entrapment in local optima. 


6.2.2 Ionic clusters 

Oleksy et al. iMi calculated the equilibrium geometries and dissociation ener¬ 
gies of electronic ground-state of He^ clusters (N = 3-35) via an extended GA 
method and employing a semiempirical valence-bond model of intracluster interac¬ 
tions. They found that for the cluster sizes of N = 3, 4, 7-9, and 14-35, the positive 
charge of He^ delocalizes over a trimer ionic core, which is more or less linear and 
centrosymmetric. Eor N = 2, 10-13, a dimer ionic core develops, while for N = 5 the 
positive charge is delocalized over the whole cluster; for N = 6, a tetramer ionic core 
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is formed. Atoms outside the ionic core are neutral and are distributed among several 
solvation rings. The electronic and geometrical structures of sodium cluster anions 
11681 [Na^, n=20-57] and cationic Na+ clusters 11691 were determined by apply¬ 
ing GA - density functional theory (DFT). Structures of Indium oxide nanoclus¬ 
ters and zirconia nanoclusters have been predicted by a method combining a robust 
evolutionary algorithm with classical interatomic potential and quantum chemical 
models ifTT^fTTn . Kim et al. ina used GA to design Nanoporous Ti 02 for Low- 
Temperature Processable Photoanodes of Dye-Sensitized Solar Cells. Darwinian 
and Lamarckian schemes within evolutionary algorithms have been compared in 
the context of structure prediction of the titania (Ti 02 ) polymorphs. Lamarckism 
in natural evolution regards the effects of “inheritance of acquired characters” as 
the motive force of evolution, while the Darwinism claims that evolution is nothing 
but the cumulative processes of natural selection with random mutation and denies 
the possibility of inheritance of acquired characters. Although the mainstream of 
today’s evolutionary theory follows Darwinism, it has been found that the Lamarck¬ 
ian scheme is more successful and efficient at generating the target structures II173L 
A combination of Buckingham and LJ potential functions describe the interactions 
between different atomic species in the system 

V{rij) = + % - % (47) 

The parameters A,y, Z?,y, C,,-, and pjj are dependent on the species involved in the in¬ 
teraction. To perform a global geometry optimization of clusters resulting from the 
microsolvation of alkali metal ions (i.e., Na+, K+, and Cs+) with benzene molecules 
Marques et al. 017411 used modified LJ function for nonelectrostatic contributions. 
Zhu et al. 11751 explored all the possible stoichiometries for Mg-O system at pres¬ 
sures up to 850 GPa. They found that two extraordinary compounds Mg 02 and 
Mg 302 became thermodynamically stable at 116 GPa and 500 GPa, respectively. 
They predicted the existence of thermodynamically stable Xe-O compounds 11761 
at high pressures (XeO, Xe 02 and Xe 03 became stable at pressures above 83, 102 
and 114 GPa, respectively). Woodley et al. 11771 predicted a greater stability of 
tetrahedral and trigonal coordinations compared to the tetragonal one for zinc ox¬ 
ide clusters [(ZnO)„, n=l-32] using an evolutionary algorithm with polarizable shell 
interatomic potentials. Wang et al. 11781 successfully applied PSO for the predic¬ 
tion of elemental (Li, C, Mg, Si), binary (Si 02 , SiC, ZnO, TiH 2 , TiB 2 , M 0 B 2 ), and 
ternary (MgSi 03 and CaC 03 ) compounds in various chemical-bonding environ¬ 
ments (metallic, ionic, and covalent bonding). Li et al. 11791 obtained the geomet¬ 
rical structures of (Al 203 )„ (n = 1-7) clusters via GA coupled with the DFT. Aval- 
troni et al. Ill801 identihed peculiar low energy isomers of small clusters (C 2 AI 4 and 
CBg^) using both standard random search and GA procedures. Pal et al. 0181111821 
determined the lowest energy structure of ZnS quantum dots of different sizes by us¬ 
ing a search based on GA coupled with the density-functional tight-binding method 
(DFTB) and found a new ring-like conhgurations of ZnS quantum dots which had 
higher HOMO-LUMO gaps compared to other ZnS quantum dot structures. Pereira 
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et al. developed an evolutionary algorithm for the global minimum energy search 
for water clusters up to (H 2 O) 20 , benzene clusters up to (C 6 H 6 ) 3 o and (CgHg)^ with 
n=2-20 II183I . Do et al. 11841 identihed the energy minima of water, methanol, wa¬ 
ter + methanol, protonated water, and protonated water + methanol clusters with 
DFT combined with basin hopping. Addicoat et al. II185II optimized a GA in or¬ 
der to identify the minimum energy structure of arbitrarily hydrogenated and hy- 
droxylated fullerenes, and implemented the method for exhaustive calculations on 
all possible isomers. They suggested that crossover operator did not have signifi¬ 
cant effect on the search efficiency of GA and hence the most efficient version of 
their GA does not employ crossover, thereby reducing it to an Evolutionary Al¬ 
gorithm (EA). Cai et al. Ml861 used GA to optimize structures of silicon clusters. 
Using the GA, Zhang et al. Ml871 predicted the most stable structures and a number 
of low-energy metastable structures for Si [001] symmetric tilted grain boundaries 
with various tilt angles which are found to be in very good agreement with the re¬ 
sults of hrst-principles calculations. Yao et al. Ml881 investigated the high-pressure 
structures of solid nitrogen through GA combined with hrst-principles electronic 
structure calculations. Hooper et al. M189I studied generic defect association com¬ 
plexes in metal oxide materials at experimentally relevant dopant concentrations by 
a specialized GA-inspired search procedure for doped metal oxides. The search al¬ 
gorithms had been tested on lanthanide-doped ceria (L = Sm, Gd, Lu) with various 
dopant concentrations. 


6.2.3 Metal Clusters & Nanoalloys 

Chen et al. M190I proposed a parallel differential evolution for cluster optimization 
(PDECO) with triangle mutation and migration operators and applied it to Pt clus¬ 
ters with great efficiency. Rogan et al. M191I implemented a GA based methody on 
small Pd Clusters. The lowest-lying isomers of the copper cluster, Cug, have been 
obtained by combining a GA driven approach with DFT M192I . Assadollahzadeh et 
al. M193I employed a seeded GA technique using DFT together with a relativistic 
pseudo-potential to search for global and energetically low-lying minimum energy 
structures of neutral gold clusters Au„ (n = 2-20). GA coupled with a tight-binding 
potential is employed by Ping et al. M194l to optimize neutral lead clusters Pb„ (n = 
2-20). Pereira et al. M195I made a study on the effectiveness of different crossover 
operators in the global optimization of atomic clusters. They came up with a ‘Cut 
and Splice’ crossover operator which proved to be extremely useful. The hrst parent 
string is “Cut” randomly along a random horizontal cutting plane into two comple¬ 
mentary parts and the second parent string is “Cut” in such a way that the number of 
the atoms beneath and above the cutting plane are equal to the number of atoms of 
the two complementary parts of the first parent. “Splice” joins the head of one string 
with the tail of the other one. Pereira et al. lisa modified this operator by the way 
it determines the sub-clusters to be exchanged. Here a random atom from the parent 
string is selected and placed on the offspring. A random number (M G — 2]) 
is generated (N is the number of atoms in the cluster). M atoms from the hrst par- 
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Fig. 5 To make the operation consistent and produce descendants with same number of atoms: the 
center of mass of the parent structures have to be translated to the origin of the coordinate system 
and have to rotate the plane of the parent clusters until the generated offspring contains the correct 
number of atoms. 


ent string closer to the first selected atom are chosen and placed on the offspring. 
The remaining (N-M-1) atoms are taken from another parent string closer to the 
first selected atom. Those atoms which are already occupied are skipped. Figure 
[^schematically depicts ‘Cut and splice crossover’. The energetic ground states of 
gold clusters, for some magic numbers, with up to 318 atoms were obtained by 
minima hopping method 11961 . A GA-DFT was invoked by Sai et al. II 1971 who 
explored the size evolution of structural and electronic properties of neutral gallium 
clusters of 20-40 atoms. The atomic structures and electronic states of Ga„ clusters 
significantly differ from the solid but resemble solid and liquid to a certain extent. 
Meng et al. 11981 found that NbCls and ZnMg intercalation in graphite results in p- 
and n-type doping, respectively. To model clusters they used Gupta semiempirical 
potential and obtained globally optimal structures of those metal clusters employing 
GA with “cut-and-splice” crossover and triangle mutation. The optimized clusters 
were then placed between two layers of graphite, each of which contained 12x12 
primitive cell and 288 atoms. The whole system was optimized with Universal po¬ 
tential and the electronic structures at the optimized doped cluster geometries were 
finally calculated using DFT. The atomic arrangement has an immense contribution 
to the catalytic activities of a nanoalloy. So, the development of a reliable method 
for predicting atomic arrangements has become increasingly important. Hong et al. 
11991 combined GA with first-principles electronic structure calculations to find the 
most stable configurations of elementary Au^ and Ag„ clusters, as well as binary 
clusters Au„,Ag„ (5 <m + n < 12). The search led to the critical Au:Ag ratios for 
which the 2D-3D transition takes place. The study showed that Ag atoms prefer pe¬ 
ripheral positions with lower coordination number while Au atoms tend to occupy 
central sites with higher coordination numbers. Johnston et al. 112001 investigated 
small Sn-Bi clusters employing a DFT/GA method. Oh et al. 12011 developed a 
combination of GA and classical molecular dynamics (MD) simulations to evalu¬ 
ate the optimal arrangement of multilayered core-shell structure for the Pt-Cu nano 
alloy with 1,654 atoms, regardless of the composition ratio. Using Tabu search in 
descriptor space and DFT Orel et al. 12021 found the global minima of the neutral 
binary Sn„,Pb„ atomic clusters, 7 < m-bn < 12, for all the possible stoichiometric 
ratios. Fournier et al. 12031 found the minimum-energy structures of Ag„Rb„ (n = 2- 
10) clusters by a combination of DFT and Tabu search. Chu et al. 12041 described a 
new fragment-based evolutionary algorithm (FA) for de novo optimization, specifi- 
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cally developed to handle organometallic and transition metal compounds. Au atoms 
preferentially segregates on the surface of Pd-Au nanoalloys. DFT calculations on 
fcc-type cubo-octahedral Pd-Au nanoparticles have indicated that the surface Pd 
atoms occupy (111) rather than (100) faces, thereby maximizing the number of 
relatively strong surface Pd-Au bonds 12051 . Larger cohesive energy, lower elec¬ 
tronegativity and smaller atomic radius (minimize bulk elastic strain) of Pd favors 
a Pd core and Au surface shell by lowering surface energy of the cluster. Besides 
electron transfer from Pd to more electro-negative Au, greater order of the Pd-Au 
than the Pd-Pd and Au-Au bonds favors Pd-Au mixing ||9^ . Ferrando et al. had 
considered a methodology based on extensive global-optimization within empirical 
potential models and subsequent DFT based local relaxation of low-energy struc¬ 
tures pertaining to different structural motifs (or basins on the energy landscape) 
of gas-phase alloy nanoclusters 111 17112071120811209112101121II12061 modelled by 
the many-body Gupta potential. Pittaway et al. 1961 performed global optimization 
of Pd-Au bimetallic clusters using a GA, coupled with the Gupta many-body em¬ 
pirical potential (EP) to model inter-atomic interactions. Johnston and coworkers 

ilTool fTOTI M Ell EHl EH Ell EH El El El E3 ES ESI did 

a series of development on theoretical simulation of alloys and bimetallic clus¬ 
ters (e.g., Cu-Ag, Ni-Ag, Cu-Au, Ag-Pd, Ag-Au, Au-Pd) with GA coupled with 
the Gupta many-body potential. AgcorePdsheii nano catalysts could be important in 
the development of micro polymer electrolyte membrane fuel cells for portable 
devices, and could also be applied to the promotion of other catalytic reactions 
Ag^Pdio cluster was identified as a wedge-shaped nano-shell with Cs point 
group symmetry by a combination of GA global optimization and DFT calcula¬ 
tion. Wu et al. 11091 performed Global optimization of Ag,„Pdn{m + n = 15) and 
Ag 3 mPdj,s- 3 m{fn = 1 — 12) clusters and came to the conclusion that silver atoms 
had a strong tendency towards segregating at the surface. Combined experimen¬ 
tal/theoretical studies have indicated the preference for PdcoreAusheit (~ 5nm) and an 
“onion-like” Pd — Au — Pd configuration for larger Pt/ —Am nanoalloys (~ I2nm). 
Johnston et al. have studied mixed coinage metal clusters, using the Gupta 

potential. 


6.3 Crystal Structure Prediction 

Evolutionary algorithms ISi EH El EH E3 EH ED EH ED EQI 
ED EH EH Eni EH EH El EH Ell have provided reliable and powerful 
means of exploring the PES of crystals leading to the identification of the most sta¬ 
ble, interesting and sometimes unexpected crystal structures. There are several other 
methods such as basin hopping ll229ll227ll2281fT42lll43l . metadynamics Il230ll23l1l . 
PSO 01781 1232 II . and simulated annealing 02331 for performing similar search. For 
large dimensional problems, due to the complex nature of PES and the presence of 
high-energy disordered structures (random initial population) these algorithms of¬ 
ten encounter serious challenges. Lyakhov et al. introduced an additional variation 
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operator - the coordinate mutation or soft mutation ll^ for the purpose of crystal 
structure prediction. Instead of using complete randomness in the mutation operator 
concerted mutation which directs the system to choices that have a higher probabil¬ 
ity to improve the fitness of the solution string is introduced. The idea of soft-mode 
mutation or softmutation 11^11341 is to move the atoms along the softest modes i.e., 
the lowest frequency eigenmodes that correspond to directions of lowest curvatures 
of the energy surface. The eigenvector corresponding to lowest non-zero eigenvalue 
determine the direction of softmutation. To calculate the softest modes one has to 
construct the dynamical matrix. It is often enough to have an approximate direction 
and sufficiently large mutation amplitude to arrive at a new low-energy structure. 
There is no need therefore of doing computationally expensive ab initio dynamical 
matrix calculation. Cheaper method for example dynamical matrix computed from 
bond hardness coefficients can be used. The soft mutation favor atoms with higher 
local order, Lyakhov et al. GSl constructed the first generation using pseudo¬ 
subcells with fractional atomic occupancies to improve the order and diversity of 
the structures. They used a fingerprint function that improves the selection process 
through removing clones. Superconducting high pressure phase of germane II135L 
structural characterization of compressed silane 12341 . rhombohedral structure of 
superhard BC 2 N 02351 . superconducting structures of BC 5 02360 . superhard mono¬ 
clinic polymorph of carbon 02371 . high-pressure orthorhombic polymorph of MgB 2 
stable above 190 GPa 02381 . high-pressure phases of CaLi 2 02390 . metallic struc¬ 
tures of oxygen at pressures in the range of 100-250 GPa 01370 . two new hexagonal 
ultra-hard phases of WN 2 02400 . two unique high-pressure metallic phases of Stan- 
nane (SnH 4 ) at superconducting temperatures of 15-22 K for the Ama2 phase at 120 
GPa and 52-62 K for the P 63 /mmc phase at 200 GPa 02411 and SiH 4 (H 2)2 at high 
superconducting temperatures of 98-107 K at 250 GPa 024211241 1 are some of the 
systems explored through ab initio evolutionary methodology. By combining PSO 
with first principles calculations, Xiang et al. 02430 discovered a new metastable Si 
phase (Si 2 o-T has a quasidirect gap of 1.55 eV) which could be a promising solar 
energy absorber. Oganov et al. Una found a new boron phase (comprised of icosa- 
hedral B 12 clusters and B 2 pairs in a NaCl-type arrangement), stable between 19 
and 89 GPa, and exhibiting evidence for charge transfer. In the context of reliably 
predicting ground state geometry, the Universal Structure Predictor: Evolutionary 
Xtallography (USPEX) Eli] [131 [l33l HH] |M1 method has turned out to be ex¬ 
tremely successful. Oganov et al. using USPEX methodology first demonstrated 
the possibility of computing hardness (at least for insulators and semiconductors) 
just from the crystal structure opens up the possibility of global optimization of 
hardness, aimed at the computational discovery of new superhard materials. They 
had introduced local measures of the quality of structure to locate defective regions 
in the crystal and used fingerprint niching, soft-mode mutation, symmetry- and 
pseudosymmetry-enabled generation of structures to greatly speed up the search for 
the global minimum 0245112461 . In a recent article 0247 1 they have summarized the 
different applications of the USPEX method as a tool for crystal structure prediction 
and showed that this method has an enormous applications in both computational 
materials design and studies of matter at extreme conditions. In a recent application 
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Oganov et al. showed in conjunction with the first-principles calculations, evolution¬ 
ary algorithm will lead to the discovery of novel dielectrics 12481 . In another appli¬ 
cation they successfully predicted energetically lower novel 2D boron structure than 
the a-sheet structure B249I . stable hafnium carbides II250II . Woodley et al. predicted 
ground state geometry for a series of transparent conducting indium oxide 11701 
and zirconia nanoclusters II171II by combining evolutionary algorithm with classical 
interatomic potential and quantum chemical models. Wu et al. showed that structure 
exploration by classical potentials with the accuracy of density functional theory is 
an efficient scheme for complex crystal structure prediction 125 1 1 . In another work 
they favorably compared Lamarckian schemes within evolutionary algorithms with 
Darwinian schemes in the context of minimum energy structure prediction of titania 
phases II73I . O’Keeffe 12521 describes several difficulties encounted in the crystal 
structure prediction. Study on GA and minima hopping method based crystal struc¬ 
ture prediction reveals that for relatively smaller sizes both methods show compa¬ 
rable efficiency while for larger systems GA becomes advantageous over minima 
hopping 022911 . Oganov et al. explained how and why evolutionary crystal struc¬ 
ture prediction works the way they do 12531 . Using adaptive GA Zhao et al. 12541 
studied the structures and stabilities of the alkaline earth metal peroxide XO 2 . They 
also predicted complex crystal structures of the orthorhombic, rhombohedral, and 
hexagonal polymorphs close to the Zr 2 Coii intermetallic compound 02551 . Li et al. 
02561 explored the high-pressure crystal structures of Mg by the PSO algorithm. 
Wang et al. have developed a software package titled Crystal Structure Analysis 
by Particle Swarm Optimization (CALYPSO), enforcing symmetry constraints on 
structure generation, bond characterization matrix for elimination of similar struc¬ 
tures, introducing partial random structures per generation for enhancing structural 
diversity, and penalty function, etc. 0178112321 for predicting crystal structure from 
random initial starting geometries. They have applied CALYPSO to crystal structure 
prediction, earth and planetary materials 0257112581125911260112611126211263112640 . 
Luo et al. 02650 predicted new stable structures of 2D boron-carbon compounds for 
a wide range of boron concentrations by the PSO algorithm implemented in the 
CALYPSO code. XtalOpt 02660 and EVO 02231 are the two other packages that 
have proved to be useful in this context. Zhu et al. proposed a method based on 
metadynamics and evolutionary algorithms 02301 and found stable and metastable 
states for Al 2 Si 05 , Si 02 , MgSi 03 , and carbon clusters 02310 from reasonable ini¬ 
tial structures providing insight into the mechanisms of phase transitions. They de¬ 
signed an evolutionary algorithm based method 02670 to automatically explore low 
energy surface reconstructions with variable surface atoms and reconstruction cells 
and illustrated it by the identification of N 3 trimeric reconstruction of GaN(lOil) 
surfaces. Liu 02681 constructed a Multi-algorithm-collaborative Universal Structure- 
prediction Environment (Muse) to efficiently find the stable and metastable struc¬ 
tures of materials under given conditions. In Muse the evolutionary algorithm was 
coupled with the simulated annealing and the basin hopping algorithms. With the 
inclusion of two new variation operators, slip and twist the performance of Muse 
was greatly enhanced. Ladda et al. 02691 described an evolutionary algorithm based 
on symmetry-preserving and symmetry-breaking mutations for the exploration of 
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the space of the conjugacy classes of crystal structures. Johnston et al. have imple¬ 
mented and parameterized both Darwinian and Lamarkian GA search to identify 
structure of clusters from experimental STEM images II270L Meredig et al. EzD 
had given considerable theoretical and computational efforts for the development of 
hrst-principles-assisted structure solution (EPASS) that automatically solve crystal 
structures. Nguyen et al. II272I performed genetic algorithm to hnd two meta-stable 
Si-IX phase with good experimental agreement in the lattice parameters. 


7 Future Directions 

The evolutionary computing techniques and their hybrids have registered impres¬ 
sive success in the elucidation of minimum energy structures of atomic clusters 
and crystalline solids. These methods have not benn explored that extensively for 
weakly bonded molecular clusters and crystals. We anticipate rapid growth of hy¬ 
brid EC based exploration of such structures. Since these surfaces are dominated by 
shallow minima separated by small barriers, special mutation and crossover opera¬ 
tors would have to be designed for exhaustive exploration of the PES and locating 
global as well as local close-lying minima. The information on such structures and 
their energies could then lead to the correct prediction of thermally averaged prop¬ 
erties of such species. 

The use of EC for directly solving the molecular Schrodinger equation has been 
very limited. The density based method described in reference ll56l has the potential 
to be a viable, even a superior alternative for computing the equilibrium molecular 
structure and the corresponding one electron density resolved in a given basis. We 
anticipate extensive exploration along that line in the future years. Intelligent inter¬ 
facing of the standard electronic structure codes and EC codes could prove instru¬ 
mental for further growth of EC techniques in exploration of molecular structures. 
One area of great promise where EC could decisively prove superior is in designing 
molecules and materials with targetted properties. Such problems can be cast in the 
mould of multimodal optimization problem which can be efficiently handled by 

the EC techniques using for example, the prey-predator model for multi-objective 
constrained optimization. Alongwith artihcial neural network, EA can meet tate-of- 
the-art methods for powerful quantitative structure-property relationships modeling 
0273112741 . We anticipate rapid growth in such studies in the next decade. 


8 Conclusions 

The last decade has seen impressive growth of EC techniques for exploring PES 
of clusters and crystals. The review highlights the important role played by soft- 
computing techniques, especially the EC methods in the elucidation of large scale 
structures on complex potential energy landscapes. When combined with DET, such 
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methods can often produce surprising results. The evolutionary computing tech¬ 
niques or for that matter the whole gamut of Soft-computing methods have not yet 
been explored exhaustively in the context of solving Schrodinger equation for even 
few electron systems although their potentials as viable tools have been demon¬ 
strated. One anticipates that the cross-talk between EC and computational chemistry 
would be productive and new hybrid algorithms for exploring electronic structures 
of atoms, molecules and clusters will evolve out of it. 
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